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ABSTRACT 


In  this  thesis  several  incompressible  oscillatory  flow  and  flutter  problems 
were  investigated.  A  previously  developed  unsteady  panel  code  for  single  airfoil 
bending  torsion  flutter  analysis  was  compared  to  Theodorsen's  classical  theory. 
The  panel  code  agrees  with  Theodorsen's  bending-torsion  flutter  analysis  for 
natural  frequency  ratios  (o)h/coa)  less  than  1.2.  In  addition,  a  two  airfoil  unsteady 

panel  code  was  modified  for  one  degree  of  freedom  flutter  analysis.  Code 
verification  was  accomplished  by  first  comparing  flat  plate  theory  to  the  unsteady 
aerodynamic  force  and  moment  coefficients  and  then  using  the  equation  of 
motion  to  determine  regions  of  instability.  The  possibility  of  active  flutter  control 
was  investigated  by  positioning  a  small  control  airfoil  in  front  of  a  neutrally  stable 
reference  airfoil.  Results  show  that  the  flutter  boundary  may  be  changed  through 
the  scaling,  placement,  or  oscillation  of  a  second  airfoil  upstream.  A  comparison 
with  pitch  damping  curves  published  by  Loewy  confirms  that  the  code  is  capable 
of  predicting  wake-induced  airfoil  flutter. 
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I.  INTRODUCTION 


A.  GENERAL 

In  this  thesis,  two  unsteady  panel  codes  UPOTFLUT  described  by  Riester 
[Ref.  1]  and  USPOTF2C  developed  by  Pang  [Ref.  2]  are  modified  and  verified. 

The  UPOTFLUT  (Unsteady  Potential  Flow  and  Flutter)  was  developed  by 
Teng  [Ref.  3]  for  unsteady  inviscid  and  incompressible  flow  over  a  single  airfoil. 
The  code  is  based  upon  the  panel  method  by  Hess  and  Smith  [Ref.  4]  that 
includes  simple  harmonic  motion  of  the  airfoil  which  is  continuously  shedding 
vortices  into  the  wake.  Riester  [Ref.  1]  extended  the  single  airfoil  code  for 
bending-torsion  flutter  analysis. 

The  USPOTF2F  (Unsteady  Potential  Flow)  code  [Ref.  2]  is  a  two  airfoil  code 
that  permits  steady  and  unsteady  analysis  of  airfoil-wake  interaction  between 
two  airfoils.  A  phase  subroutine  is  added  to  convert  the  time  dependent  lift  and 
moment  histories  into  the  frequency  domain. 

Building  upon  the  work  of  Riester  [Ref.  1]  which  focused  mainly  on 
aerodynamic  verification  of  UPOTFLUT  code,  the  two  dimensional  bending- 
torsion  flutter  problem  is  examined.  Theodorsen  [Ref.  5,  page  11]  presents  flat 
plate  flutter  speeds  for  comparison  to  UPOTFLUT  computations.  Specific  cases 
are  selected  to  compare  the  aerodynamic  forces  and  solutions  of  the  flutter 
determinant. 

The  two  airfoil  USPOTF2  code’s  unsteady  aerodynamics  are  verified  for  pitch 
and  plunge.  Then  the  instability  of  one  degree  of  freedom  pitching  oscillations  is 
chosen  to  explore  active  flutter  control.  Pitch  damping  dependence  upon 
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upstream  airfoil  size  and  position  is  fully  explored.  Emphasis  is  placed  upon  the 
wake  effect  and  the  code's  ability  to  predict  wake-induced  flutter. 

Whenever  possible,  the  theoretical  work  of  Theodorsen  [Ref.  5],  Smilg  [Ref. 
6],  Loewy  [Ref.  7]  and  other  numerical  examples  are  used  for  comparison  to 
verify  accuracy  and  trends.  However,  for  many  of  the  computations  performed 
no  other  analytical  theories  and  applicable  experimental  data  exist. 

B.  SCOPE 

Chapter  II  contains  bending-torsion  flutter  verification  using  the  single  airfoil 
UPOTFLUT  code  and  Theodorsen  theory.  Aerodynamic  verification  of  force  and 
moment  coefficients  using  the  two  airfoil  USPOTF2F  code  is  accomplished  in 
Chapter  III.  Chapter  IV  develops  the  equations  of  motion  for  steady-state 
oscillations  in  one  degree  of  freedom.  An  in-ground  effect  numerical  simulation  is 
completed  by  oscillating  airfoils  out  of  phase.  The  positioning  and  size  of  the 
control  airfoil  for  active  flutter  control  are  investigated  in  Chapter  V.  Application 
to  rotary  wings  along  with  a  comparison  to  wake-induced  flutter  theory  is  also 
discussed. 
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II.  SINGLE  AIRFOIL  ANALYSIS 


The  results  discussed  and  presented  in  this  chapter  were  produced  from  the 
UPOTFLUT  code.  This  single  airfoil  analysis  builds  on  work  by  Riester  [Ref.  1]. 

A.  INTRODUCTION  TO  UNSTEADY  PANEL  CODE  THEORY 

Given  a  two-dimensional  airfoil,  the  governing  equation  to  solve  the  inviscid, 
incompressible  flow  over  the  airfoil  is  Laplace’s  equation.  Since  the  two- 
dimensional  Laplace  equation  is  linear,  the  principle  of  superposition  may  be 
applied.  The  four  elementary  flows:  uniform  flow,  source  flow,  doublet  flow  and 
vortex  flow  are  used  to  determine  the  flow  around  arbitrary  bodies.  For  non¬ 
lifting  bodies  a  uniform  flow  and  source  flow  are  required.  For  lifting  bodies  the 
vortex  flow,  must  also  be  included  to  provide  circulation.  The  boundary 
conditions  which  must  be  satisfied  are  flow  tangency  at  the  surface  and  the  Kutta 
condition  at  the  trailing  edge. 

Steady-state  panel  methods  as  described  by  Hess  and  Smith  [Ref.  4]  may  be 
extended  to  airfoils  in  simple  harmonic  motion.  The  continuous  change  in  the  lift 
implies  a  continually  changing  circulation  about  the  airfoil.  The  Helmholtz  vortex 
theorem  requires  that  any  change  in  the  circulation  around  an  airfoil  must  be 
matched  by  the  appearance  of  an  equal  counter-vortex  or  starting  vortex  to 
achieve  constant  total  circulation  in  the  flow  field.  The  starting  vortices  are 
assumed  to  be  shed  from  the  trailing  edge  after  each  time  step  and  move  with  the 
local  flow  velocity. 

Riester  [Ref.  1],  page  23,  showed  that  for  single  airfoil  analysis  3  cycles  is  a 
sufficient  wake  length  for  amplitude  and  phase  analysis  when  Kp  =  1 .0.  The 
benchmark  chosen  for  computations  is  3  cycles,  65  time  steps  per  cycle.  A 
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NACA0007  airfoil  geometry  was  chosen  for  comparison  to  flat  plate  theory.  A 
thinner  airfoil  does  not  provide  accurate  results  due  to  insufficient  number  of 
panels  to  adequately  resolve  the  leading  edge  and  vortex  interaction  between 
top  and  bottom  surfaces. 

B.  TWO  DEGREE  OF  FREEDOM  EQUATIONS  OF  MOTION  AND 

FLUTTER  DETERMINANT 

If  a  positive  concentrated  load,  on  a  cantilever  wing  of  uniform  cross  section 
is  applied  at  the  leading  edge  of  the  wing  tip,  the  tip  cross  section  will  be 
displaced  upwards  and  twist,  thus  increasing  the  angle  of  attack.  If  a  similar 
concentrated  load  is  applied  at  the  trailing  edge,  the  tip  cross  section  will  flex 
upwards  and  the  angle  of  attack  will  decrease.  A  point  exists  between  the 
leading  edge  and  trailing  edge  where  a  concentrated  load  may  be  applied  without 
twisting  or  causing  the  airfoil  to  rotate.  This  point  is  called  the  shear  center  or  the 
flexural  center  of  the  airfoil.  The  elastic  axis  is  defined  as  the  locus  of  shear 
centers  of  the  cross  section,  thus  it  is  a  natural  reference  line. 

Forces  acting  upon  the  airfoil  are  the  elastic  restoring  forces,  inertial  forces 
and  aerodynamic  forces.  Displacement  variables  for  two  degree  of  freedom, 
simple  harmonic  motion  are  pitch  (a)  and  plunge  (h): 

Ca  2  torsional  stiffnep.!  of  wing 
Ch  2  stiffness  of  wing  in  flexure 
M  2  mass  of  wing  per  unit  span 
Sa  2  static  moment  referred  to  elastic  axis 
Ia  2  moment  of  inertia  about  elastic  axis 
The  elastic  restoring  forces  are  hCj,  and  aCa.  The  inertial  forces  respectively  in 
pitch  and  plunge  are: 

a  Ia  +  h  Sa 
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h  M  +  a  Sa 

Applying  equilibrium  by  equating  the  inertial  forces  to  the  external  aerodynamic 
and  structural  elastic  forces  and  assuming  zero  structural  damping  yields: 

a  Ia  +  h  Sa  =  Moment  -  aCa 
h  M  +  a  Sa  =  Lift  -  hCh 
The  natural  frequencies  in  pitch  and  plunge  are: 

(0h  = 

(0o  = 

Substituting  into  the  equations  of  motion: 

hM  +  aS„  +  hMtoJ  =  Lift 
ala  +  h  Sa  +  aIo0)£  =  Moment 

For  simple  harmonic  motion: 

h  =  h0eitot  a  =  a0eia>t 
h  =  -a>2hoe,t0t  a  =  -a>2aoe1{0t 
e“"*  (-to2  h0  M  -  co2  a0  Sa  +  h0  McoJ )  =  Lift 
ei<at  (-co2a0Ia  -co2h0Sa  +  a0Ia£o2 )  =  Moment 

The  panel  code  nondimensionalizes  lift  and  moment  using  the  characteristic 
length  of  chord  (2b)  vice  semichord. 

Lift  =  (CLo  +Cu,)ipU22b 

Moment  =  (CMa  +  CMh )  |  pU2  4  b2 

Introducing  dimensionless  flutter  parameters: 

Wing  density  parameter  u.  =  ■  —5- 

Ttpb 
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Center  of  gravity  from  elastic  axis:  xa  =  -~ 

Mb 

Radius  of  gyration:  r*  =  — 

Mb 


Flutter  frequency  ratio:  X = j 
The  equations  of  motion  become: 

^  [(-to*  h„  M  -  a>2  a„  Sa  +  h0  Mo>2 )  =  (C^  +  Cu, )  ±  pU2  2  b] 
^2  [(-<o2«„I„  -  to2  h0  Sa  +  a0Ioa>|)  =  (CMa  +Cm  )ipU24b2j 
Jk_„  x  +M^Y=(C  +C  )WU-**L 

b  oa  blcoj  1x1  121  itMco2  4b2 

r*  - £»  x  +  a  r2  =  (C  +C  1 

*  °  b  “  °  a  l  (0  J  (  m“  “  J  7tM(024b2 

-a0x0+|-^)  -l]-(Cu+Cu)^ 

“<>r“((lf)‘1)"'b'^o)  =  (CMa  +C“)^^ 

Rearrange  the  two  equations  into  the  form: 

A®  +  B(a)  =  0 

D®  +  E(a)  =  0 


where 


l  UJUJJ  i) 


B  =  11X0  +  4^2 

jcaKp 
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D  =  (ixa  + 


4C 


Mh 


E  =  |irj  1-  -*■ 


TtaKp 


A  B 
|D  E 


+  A 


lm  aginary 


To  find  the  critical  flutter  speed  and  frequency  the  characteristic  equation  is 
solved  for  U  and  to.  Since  A  is  complex  both  the  real  and  imaginary  parts  must 
equal  zero  yielding  two  real  quadratic  equations  for  two  unknowns.  The  real 
equation  typically  assumes  a  parabolic  shape  and  the  imaginary  equation 
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C.  MODIFICATION  OF  UPOTFLUT 

Significant  modifications  to  the  UPOTFLUT  code  involved 
nondimensionalizing  flutter  determinant  parameters. 

The  phase  subroutine  converts  the  time  history  of  lift  and  moment  into  the 
form: 

F(t)=Amp*  cos(a>t+<(>) 

Where  Amp  is  the  amplitude,  (0  is  frequency  of  oscillation  and  <j>  is  igle  of 
phase  lag  of  the  aerodynamic  forces  to  the  motion  of  the  airfoil.  An  iterative 
numerical  scheme  to  find  the  phase  shift  involves  summing  the  differences 
between  actual  time  history  and  the  fitted  curve  over  a  half  cycle,  then  shifting  <)> 

to  minimize  error.  When  iterations  start  with  a  phase  error  of  n  the  total  error  is  a 
maximum,  hence  the  error  curve’s  derivative  equals  zero.  If  initial  phase  shift  step 
size  is  too  small,  the  numerical  scheme  will  remain  at  a  maximum  instead  of 
converging  to  minimize  error.  Initial  step  size  of  four  degrees  is  sufficient  to 
ensure  convergence. 

D.  VERIFICATION  OF  UPOTFLUT 

NACA  TR-685  page  1 1  [Ref.  5]  provides  the  closed  form  solution  of  a  flat 
plate  for  comparison  with  the  UPOTFLUT  code.  The  results  of  Case  (a)  are 
identical  to  code  validation  performed  by  Riester  [Ref.  1]  page  101.  Case  (a)  and 
Case  (b)  results  are  tabulated  in  Table  2.1.  In  general,  the  panel  code  results  agree 
with  Theodorsen’s  theory.  Next,  the  xalpha  =  0.2  Case  (h)  was  calculated  using 
UPOTFLUT  panel  code.  Case  (h)  physically  represents  a  heavy  wing  with  the 
center  of  mass  aft  of  the  elastic  axis.  The  pivot  axis  is  located  at  35  percent  of 
chord.  These  parameters  make  the  cross  section  inherently  susceptible  to  flutter. 
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As  illustrated  in  Figure  2.3,  the  panel  code  results  and  Theodorsen  theory  diverge 
at  natural  frequency  ratios  greater  than  0.8. 

The  Theodorsen  theoretical  curve  was  verified  using  two  dimensional 
flutter  analysis  described  by  Fung,  page  235  [Ref.  8].  The  analysis  was 
completed  by  programming  the  fundamentals  of  flutter  theory  in  MATLAB 
[Appendix  A]  using  unsteady  aerodynamic  coefficients  found  on  page  412  of 
Scanlan  and  Rosenbaum  [Ref.  9].  The  non-linear  complex  quadratic  equations  are 
solved  by  using  a  bisection  root  solving  numerical  scheme.  Calculated  non- 
dimensional  flutter  speeds  coincided  with  published  coefficients  for  Q>h/(oa  less 
than  1.2.  For  coh/<o0  greater  than  1.2,  the  critical  reduced  frequencies 
encountered  are  greater  than  4.0  based  upon  semichord  (Kt  =4.0).  At  high 

reduced  frequencies,  the  flutter  determinant  becomes  ill  conditioned  and  very 
sensitive  to  small  changes  in  any  of  the  flutter  variables.  Theodorsen  did  not 
have  the  computational  power  of  computers  in  1940  when  TR-685  was 
published.  Therefore,  his  results  were  recomputed.  Although  Theodorsen’s 
curve  for  Case  (h)  in  TR-685,  page  1 1  was  found  to  be  inaccurate  for  frequency 
ratios  greater  than  1.25  (as  shown  in  Figure  2.3),  the  new  theoretical  curve  using 
numerical  methods  did  not  conect  the  large  discrepancy  between  his  theory  and 
the  panel  code.  The  source  for  this  discrepancy  may  come  from  the  computation 
of  the  unsteady  aerodynamics  or  from  the  flutter  calculations.  Figures  2.4 
through  2.11  and  Table  2.3  indicate  close  agreement  between  the  aerodynamic 
coefficients  for  lower  reduced  frequencies  with  errors  increasing  slightly  at  higher 
frequencies.  The  panel  code  unsteady  aerodynamic  coefficients  are  derived  and 
discussed  by  Riester  [Ref.  1]  and  reproduced  here  specifically  for  TR-685  page  1 1 
Case  (h). 
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The  agreement  between  Theodorsen’s  and  the  panel  code  unsteady 

aerodynamics  suggests  that  the  discrepancy  in  flutter  coefficient  is  mainly  due  to 

the  solution  of  the  flutter  determinant.  The  flutter  coefficient  is  defined  as: 

U  1  1 

bo)„  ”  VX  K, 

Figure’s  2.12  through  2.14  and  Table  2.4  compare  the  real  and  imaginary 
roots  of  the  flutter  determinant  from  theory  and  the  panel  code.  The  intersection 
point  of  the  real  and  imaginary  roots  represents  the  flutter  frequency.  It  is  shown 
in  Figure  2.13  that  the  critical  reduced  frequency  may  differ  by  up  to  a  factor  of 
three  from  Theodorsen’s  result  for  high  critical  frequencies  and  high  frequency 
ratios. 
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AIRFOIL  TYPE  s  NACA  0007  AIRFOIL 
N LONER  -  100  ,  N UPPER  •  100 


I Ft AO  N LOWER  NUPPER 
0  100  100 
AIRFOIL  TYPE 
3 

IRAHP  IOSCIL  ALPI  ALPMAX  PIVOT 

0  1  -1.0  1.0  0.0 

FREQ  RFOSTP  RFQFNL 
0.04  0.02  0.2 

IGUST  UGUST  VGUST 
0  0.  0. 

ITRANS  DELHX  OELHY  DELI  PHASE 
0  0.00  .01  -.01  0.00 
CYCLE  NTCYCLE  TOL 

3  65  0.001 

naot  t  naot  X  aoa  values  Multiplied  by  10  (integer) 

2  05  10  20  25  39  50 

wh/walpha  ratio  Maes  ratio  xalpha  (ralpha)**2 

0.4  3  0.2  0.25 

ComMenta . . . 

I RAMP  0:  n/a  RFREQ  is  based  on  full  chord 

1:  straight  rasp 
2:  Modified  rasp 

IOSCIL  Ot  n/a  RFREQ  is  based  on  full  chord 

1:  Sinusoidal  pitch.  Motion  starts  at  Min  Aoa 

pivot  :  Position  of  elastic  axis  if  leading  edge-0  and  trailing  edge-1 

ITRANS  0:  n/a 

1:  Translational  harMonic  oscillation 

CYCLE  :  /  of  cycles  for  oscillatory  Motions 

-In  case  of  rasp,  cycle-1.5  denotes  airfoil  is  held 
at  Max  aoa  for  the  duration  of  .5  cycle 
-For  steady  state  solution  set  it  to  0 

ntcycle:  #  of  tine  steps  for  each  cycle 

CYCLE*NTCYCLE  is  United  to  200  currently. 

NAOT:  #  of  input  aoa  for  co  output 

-  angles  should  be  in  Increasing  order, 

-  for  oscilatory  notions  angles  should  increase 
first,  then  decrease.  Decreasing  angles  are  for 
the  return  cycle.. 


Figure  2.2  Typical  UPOTFLUT.IN  input  file 
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TABLE  2.1  TR-685  CASE(A)  and  CASE(B)  COMPARISON 
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Flutter  Coefficient  Comparison 
(TR-685  pg.  11  Case  h  Xalpha=0.2) 


Rgure  2.3  Rutter  Coefficient  Comparison 
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Real  (Cl 


Re(Cl)  Due  to  Pitch  Comparison 
(TR  -  685,  Case  h,  pg.  11;  200 

panels;  3  cycles,  65  time  steps 
per  cycle) 


Kp  =  2bw/U 


Figure  2.4  Real(Cl)  Due  to  Pitch 
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Imaginary  (Cl) 


Imaginary  Part  of  Cl  Pitch 
Comparison  (TR  -685,  Case  h,  pg. 
11;  200  panels;  3  cycles  65  steps 

per  cycle) 


Kp  =  2bw/U 


Figure  2.5  Imaginary(Cl)  Due  to  Pitch 
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Figure  2.6  Real(Cm)  Due  to  Pitch 


Im(Cm)  Due  to  Pitch  Comparison 
(TR  -685,  Case  h,  pg.  11;  200 
panels;  3  cycles,  65  time  steps 
per  cycle) 


■o.oi 
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Re(CI)  Due  to  Plunge  Comparison 
(TR  -  685,  Case  h,  pg.  11;  200 
panels;  3  cycles,  65  time  steps 
per  cycle) 


Figure  2.8  Real(Cl)  Due  to  Plunge 
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Figure  2.9  Imaginary(Cm)  Due  to  Plunge 
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Figure  2.1 1  Imaginary(Cm)  Due  to  Plunge 


0.0822  I  -0.0688  I  0.9379  I  -0.2900  I  0.1340  I  0.0217 


SQRT(X) 


Rea!  and  Imaginary  Roots  of 
SQRT(X)  for  wh/walpha  =  0.4 


Figure  2.12  Roots  for  Frequency  Ratio  of  0.4 


Real  and  Imaginary  Roots  of 
SQRT(X)  for  wh/walpha  =  1.0 


Figure  2.13  Roots  for  Frequency  Ratio  of  1 .0 
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SQRT(X) 


Real  and  Imaginary  Roots  of 
SQRT(X)  for  wh/walpha  =  1.6 


Figure  2.14  Roots  for  Frequency  Ratio  of  1.6 
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TABLE  2.4  SOLUTIONS  OF  FLUTTER  DETERMINANT 


The  roots  SQRT(X)  of  the  real  and  imaginary  equations  against  1/Kt. 


TR  •  685  pg.  1 1  Case  h;  mass  ratio  »  4;  a  -  -0.3;  Xalpha  -  0.2;NACA0007 


itch:  0.01  h/2b  plunge;  200  panels;  3  cycles  65  time  steps 


wratio  0.4 


anel 


Reduced  Freque 


1/Ktheodorsen  Real  Panel 


0.25 


0.50 


0.83 


1.25 


1.67 


2.00 


2.50 


3.33 


5.00 


0.9628 


0.9771 


1.0273 


1.1422 


1.3145 


1.4963 


1.9345 


Imagina 


1.3573 


1.3533 


1.3386 


1.3315 


1.3338 


1.3400 


1.3536 


1.3689 


1.4982 


Real  Theo 


0.9783 


1.0047 


1.0675 


1.1899 


1.3621 


1.5403 


1.9409 


Imagina 


1.3793 


1.3723 


1.3747 


1.3820 


1.3925 


1.4025 


1.4187 


1.4465 


1.4974 


wratio  1.0 


Reduced  Freque 


1/Ktheodorsen  Real  Panel 


0.25 


.50 


0.83 


1.25 


1.67 


0.8235 


0.8293 


0.8577 


0.9223 


1.0026 


1.0602 


1.1184 


1.1628 


1.2075 


Imagina 


0.8843 


0.8743 


0.8642 


0.8607 


.8631 


.8672 


0.8752 


0.8820 


0.9556 


theor 


Real  Theory  Imagina 


0.8317 


0.8489 


0.8880 


0.9561 


1.0303 


1.0806 


1.1293 


1.1624 


1.1767 


0.8679 


0.8661 


0.8702 


0.8772 


.8852 


0.8917 


0.9013 


0.9160 


0.9394 


Reduced  Fr 


1/Ktheodorsen 


0.25 


0.50 


•  0.83 


1.25 


1.67 


2.00 


2.50 


3.33 


5.00 


Real 


5 


.5910 


0.5981 


0.6167 


0.6388 


0.6546 


0.6737 


0.6968 


0.7313 


wratio  1.6 


net 


Imaginary  Real  Imagina 


0.6148  0.5947  0.5958 


0.6059  0.6005  0.5951 


0.5988  0.6130  0.5987 


0.5966  0.6330  0.6041 


0.5985  0.6533  0.6098 


0.6014  0.6675  0.6143 


0.6068  0.6843  0.6208 


0.6107  0.7022  0.6302 


0.6592  I  0.7184  0.6442 
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III.  AERODYNAMIC  COEFFICIENT  VERIFICATION  OF 

TWO  AIRFOIL  CODE 

The  two  airfoil  USPOTF2F  code  was  modified  for  one  degree  of  freedom 
flutter  analysis.  This  chapter  verifies  the  aerodynamic  force  and  moment 
coefficients  by  comparison  to  flat  plate  theory. 

A.  INTRODUCTION  TO  TWO  AIRFOIL  PANEL  CODE 

The  USPOTF2F  unsteady  panel  code  is  essentially  the  single  airfoil  UPOT 
code  extended  to  two  airfoils.  Airfoil  number  one  is  treated  as  the  master  airfoil 
and  remains  at  the  origin  of  a  cartesian  coordinate  system.  Airfoil  number  two 
may  be  arbitrarily  placed  in  reference  to  the  master  airfoil,  as  long  as  the  leading 
airfoil’s  wake  does  not  impinge  upon  the  trailing  airfoil  or  its  wake.  All 
aerodynamic  and  flutter  calculations  pertain  to  the  master  airfoil. 

The  phase  subroutine  computes  the  complex  unsteady  aerodynamic 
coefficients  (refer  to  page  8),  then  the  subroutine  writes  the  coefficients  to 
pitch.in  and  plunge.in  files.  Piteh.in  contains  Kp,  RealfC^),  Imaginary 

Real(CMa)  and  Imaginaiy(CMo)  tabulated  in  columns  from  left  to  right  respectively 
for  single  degree  of  freedom  pitch.  Plunge.in  contains  the  same  information 
except  the  coefficients  are  due  to  single  degree  of  freedom  plunge  in  simple 
harmonic  motion.  The  phase  subroutine  searches  the  last  cycle  of  the  time  history 

of  lift  and  moment  for  a  maximum  and  minimum  value.  Average  lift  and  moment 
coefficients  are  calculated  from  (CLnux  -CLmin)/2.  Each  average  is  then 

subtracted  from  the  time  history  respectively  before  curve  fitting  to  a  sine 
function: 

F(t) = Amp  *  sin(cot + <j>) 
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A  function  of  this  form  assumes  a  zero  average,  or  equal  maximum  and  minimum 
magnitudes.  Because,  the  pitch  and  plunge  values  start  from  zero  at  time  equal  to 
zero,  thus  this  initial  condition  requires  a  sine  function.  The  phase  portion  uses 
the  first  k  radians  of  the  last  cycle  for  curve  fitting.  Selecting  the  positive  area  of 
the  sine  curve  for  integration  avoids  errors  near  k/2  and  3  jc/2  by  not  including 
positive  and  negative  values  in  the  integration  interval. 

The  product  of  cycles  and  time  steps  may  not  exceed  200  due  to  array  sizing. 
As  a  warning  to  the  user,  note  that  the  period  of  the  last  cycle  is  based  upon  the 
reduced  frequency  of  the  master  airfoil.  If  the  airfoils  are  oscillating  at  different 
reduced  frequencies,  the  force  and  moment  data  for  airfoil  two  in  pitch.in  and 
plunge.in  files  are  invalid. 

To  plot  the  wake  vortices  position  and  airfoil  geometry,  use  the  original 
USPOTF2C  code.  Fort.  10  and  fort.  11  contain  wake  vortices  position  for  each 
time  step.  Cutting  data  from  the  last  time  step  in  fort.  10  and  fort.ll,  and  then 
pasting  the  data  into  separate  files  will  allow  for  plotting  applications. 

B.  THEORETICAL  COEFFICIENT  ORIGIN 

The  theoretical  unsteady  aerodynamic  coefficient  values  for  a  flat  plate  were 
calculated  using  the  Tables  of  Subsonic  Incompressible  Aerodynamic  Coefficients 
from  Scanlan  and  Rosenbaum  [Ref.  9]  pg.  412.  Riester  [Ref.  1  pg.  18]  derives  the 
complex  coefficient  of  lift  and  moment  values.  The  constants  and  equations  were 
programmed  in  MATLAB  [Appendix  A]  for  comparison  to  USPOTF2F. 

C.  GEOMETRY  FOR  SINGLE  AIRFOIL  COMPARISON 

Airfoils  are  separated  by  50  chord  lengths  with  no  shift  in  the  horizontal 
direction.  The  interaction  between  the  two  airfoils  is  assumed  to  be  negligible  at 
this  distance.  Inputs  of  200  panels  for  each  airfoil,  pitch  and  plunge  amplitudes  of 
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one  degree  and  0.01  chord  are  used.  Riester  [Ref.  1]  showed  that  small 
amplitudes  yield  more  accurate  results,  especially  in  plunge  at  high  reduced 
frequencies.  As  with  UPOTFLUT  aerodynamic  verification,  3  cycles,  65  time 
steps  per  cycle  was  chosen  as  a  benchmark  wake  length. 

D.  AERODYNAMIC  VERIFICATION 

The  comparison  between  theory  and  USPOTF2F  code  is  presented  in  Figures 
3.3  through  3.18  and  Tables  2.1  through  2.4.  For  clarity,  the  computed 
aerodynamic  coefficients  are  presented  in  phasor  form  as  a  magnitude  and  phase 
angle  vice  complex  form.  Phasor  form  allows  for  easier  identification  of  errors. 

It  is  seen  that  the  panel  code  results  are  in  agreement  with  theory. 
Magnitude  errors  are  less  than  five  percent  and  phase  error  typically  ranged  from 
one  to  three  degrees.  A  portion  of  the  error  may  be  attributed  to  comparison  of 
an  airfoil  of  finite  thickness  to  a  theoretical  flat  plate. 
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NUMBER  OF  LINES  FOR  TITLE 
1 

TWO  NACA0007  OCSILLATIHG  AIRFOILS 

I  FLAG  N LOWER  NUPPER 

0  73  75 

NAIRFO,  XSHIFT  .  ¥ SHIFT,  SCALE 
2  0  -100  1.0 

HACA  AIRFOIL  TYPE, 

7 

7 

ALF1  ALP2  DALF1  DALF2  TCON1  TCON2 

0.0  0.0  0.0  0.0  0.0  0.0 

FREQ1  FRBQ2  PIVOT1  PIVOT2 

3.0  3.0  0.0  0.0 

UGUST  VGUST  DELHX1  DEL1IX2  OELilVl  DELHV2  PHASE1  PHASE2 
0.0  0.0  0.0  0.0  0.1  0.1  0.0  0.0 
TF  OTS1  DTS2  TOL  TADJ  SCL  S CM  SCAM  NGIES 
6  65  0.0  .001  0.0  0.0  0.0  0.0  1 

STEADY  OUTPUT 

false  false 

STEADY— TRUE  IF  ONLY  STEADY  SOLUTION.  FALSE  OTHERWISE. 

OUTPUT— TRUE  IF  YOU  WANT  COMPLETE  OUTPUT  TO  SCREEN. 

I FLAG  :  O  IF  AIRFOIL  IS  NACA  XXXX  OR  230XX 
1  OTHERWISE. 

H LOWER  :  NO.  OF  PANELS  USED  ON  BOTH  AIRFOIL  LOWER  SURFACES. 

NUPPER  :  NO.  OF  PANELS  USED  ON  BOTH  AIRFOIL  UPPER  SURFACES. 

NAIRFO  !  NUMBER  OF  AIRFOILS. 

XSHIFT  :  RELATIVE  X  DIST.  FROM  2  AIRFOIL  PIVOT  POSITION  WRT 
GLOBAL  COORDINATE  SYSTEM. 

YSlilFT  ;  RELATIVE  Y  DIST.  FROM  2  AIRFOIL  PIVOT  POSITION  WRT 
GLOBAL  COORDINATED  SYSTEM. 

HACA  AIRFOIL  TYPE  :  ENTER  NACA  4  OR  5  DIGIT  CODE  FOR  AIRFOILS., 

IF  NOT  A  NACA  AIRFOIL,  SUPPLY  AIRFOIL 
X(I) ,Y(I)  COORDS.  FOR  BOTH  AIRFOILS  IN 
FILE  CODE  2. 

ALF 1/2  :  INITIAL  ANGLE  OF  ATTACK  FOR  AIRFOILS  IN  DEGREES. 

DALP1/2  :  CHANGE  IN  AOA  IN  DEGREE  FOR  NON  OSCILL.  MOTION. 

MAX  AMPLITUDE  OF  AOA  IN  DEGREE  FOR  ROT.  HARMONIC  MOTION. 

TCON1/2  t  NON-DIMENSIONAL  RISE  TIME  (Vinf.t/C)  OF  AOA  FOR 
MODIFIED  RAMP  CHANCE  IN  AOA. 

FREQ1/2  s  NON  DIMENSIONAL  OSCILL.  (WC/Vinf.)  FOR  HARMONIC  MOTIONS. 
riVOTl/2:  LENGTH  FROM  LEADING  EDGE  TO  PIVOT  POINT  FOR  LOCAL  SYSTEM. 

(THE  GLOBAL  SYSTEM'S  ORGIN  IS  THE  FIRST  AIRFOILS  PIVOT  POSITION 
UGUST  :  HAG.  OF  NON-DIM.  GUST  VELOCITY  ALONG  GLOBAL  X  DIRECTION. 

VGUST  :  MAG.  OF  NON-DIM.  GUST  VELOCITY  ALONG  GLOBAL  Y  DIRECTION. 

DELI1X1/2 :  NON-DIM.  TRANSLATIONAL  CHORDWISE  AMPLITUDE. 

DELIIY1/2:  NON-DIM.  TRANSLATIONAL  TRANSVERSE  AMPLITUDE  (plunging). 
NIASE1/2!  PHASE  ANGLE  IN  DEGREE  BETWEEN  CHORDWISE  AND  TRANSVERSE 

TRANSLATIONAL  OSCILL.  WITH  THE  LATTER  REF.  TO  THAT  AIRFOIL. 

TF  :  FINAL  NON-DIM.  TIME  TO  TERMINATE  UNSTEADY  FLOW  SOLUTION. 

DTS  :  STARTING  TIME  STEP  FOR  NON-OCIILL  MOTIONS (TADJ-0) . 

NO.  OF  COMPUTAIONAL  STEPS  PER  CYCLE  FOR  HARMONIC  MOTION 
(FOR  2  FREQ  OCILL.  IT  USES  THE  LARGEST  FREQ) 

BASELINE  TIME  STEP  FOR  ALL  MOTIONS (TADJ  NOT  -0) 

0TS2  !  STARTING  NON-DIM  TIME  FOR  SECOND  AIRFOIL  MOTION  TO  BEGIN. 

(O  .0  BEGIN  MOTION  AT  THE  SAME  TIME) . 

TOL  :  TOLERANCE  CRIERION  FOR  CONVERGENCE  FOR  (Uw)k  and  (Vw)k. 

TADJ  :  FACTOR  BY  WHICH  DTS  WILL  BE  ADJUSTED. 

SCLA  :  STEADY  LIFT  COEFF.  FOR  THE  SINGLE  AIRFOIL  AT  THE  SPEC.  AOA. 

SCM  :  STEADY  MOMENT  COEFF.  FOR  THE  SINGLE  AIRFOIL. 

SGAH  :  STEADY  VORTICITY  STRENGTH  FOR  THE  SINGLE  AIRFOIL. 

NGIES  :  OPTION  TO  CHANGE  THE  UNSTEADY  KUTTA  CONDITION. 

0  EQUAL  PRESSURE  AT  THE  TRAILING  EDGE  PANELS. 


Figure  3.1  USPOTF2C  Input  Namelist  for  Wake  Analysis  (Original) 
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NUMBER  OF  LINES  FOR  TITLE 
1 

TWO  NACA0007  OCSILLATING  AIRFOILS 

IF! AG  H LOWER  NUPPER 
0  100  100 
NAIRFO ,  XSHIFT  ,  YSIIIFT,  SCALE 

2  -43.98  -25  1.0 

UACA  AIRFOIL  TYPE, 

7 

7 

FREOl  FREQ 2  DALP1  DALP2 

0.1  0.1  1.0  1.0 

FREQSTEP1  FREQSTEP2  PIVOT1  PIVOT2 

0.0  0.0  0.0  0.0 

FIHALFREQ1  FINALFREQ2  DELHY1  DELHY2 

0.05  0.05  0.01  0.01 

CYCLES  DTS1  0TS2  TOL  TADJ  SCL  SCM  SGAM  NGIES 

3  65  0.0  0.001  0.0  0.0  0.0  0.0  1 

STEADY  OUTPUT 

false  false 

STEADY— TRUE  IF  ONLY  STEADY  SOLUTION.  FALSE  OTHERWISE. 

OUTPUT— TRUE  IF  YOU  WANT  COMPLETE  OUTPUT  TO  SCREEN. 

IFiAG  O  IF  AIRFOIL  IS  NACA  XXXX  OR  230XX 
1  OTHERWISE. 

NLOWER  :  NO.  OF  FANEI.S  USED  ON  BOTH  AIRFOIL  LOWER  SURFACES. 

NUPPER  :  NO.  OF  PANELS  USED  ON  BOTH  AIRFOIL  UPPER  SURFACES. 

NAIRFO  :  NUMBER  OF  AIRFOILS. 

XSHIFT  !  RELATIVE  X  DIST.  FROM  2  AIRFOIL  PIVOT  POSITION  WRT 
GLOBAL  COORDINATE  SYSTEM. 

YSIIIFT  !  RELATIVE  Y  DIST.  FROM  2  AIRFOIL  PIVOT  POSITION  WRT 
GIOBAL  COORDINATED  SYSTEM. 

NACA  AIRFOIL  TYPE  S  ENTER  NACA  4  OR  5  DIGIT  CODE  FOR  AIRFOILS., 

IF  NOT  A  NACA  AIRFOIL,  SUPPLY  AIRFOIL 
X(l) , Y(I)  COORDS.  FOR  BOTH  AIRFOILS  IN 
FILE  CODE  2. 

DALP1/2  :  CHANGE  IN  AOA  IN  DEGREE  FOR  NON  OSCILL.  MOTION. 

MAX  AMPLITUDE  OF  AOA  IN  DEGREE  FOR  ROT.  HARMONIC  MOTION. 

FREQ1/2  :  NON  DIMENSIONAL  OSCILL.  (WC/Vlnf.)  FOR  HARMONIC  MOTIONS. 
PIVOT1/2:  LENGTH  FROM  LEADING  EDGE  TO  PIVOT  POINT  FOR  LOCAL  SYSTEM. 

(THE  GLOBAL  SYSTEM'S  ORGIN  IS  THE  FIRST  AIRFOILS  PIVOT  POSITIOF 
DELHY1/2:  NON-DIM.  TRANSLATIONAL  TRANSVERSE  AMPLITUDE  (plunging). 

DTS  :  STARTING  TIME  STEP  FOR  NON -OC I ILL  MOTIONS (TADJ-0) . 

NO.  OF  COHPUTAIONAL  STEPS  PER  CYCLE  FOR  HARMONIC  MOTION 
(FOR  2  FREQ  OCILL.  IT  USES  THE  LARGEST  FREQ) 

BASELINE  TIME  STEP  FOR  ALL  MOTIONS (TADS  NOT  -0) 

DTS2  :  STARTING  NON-DIM  TIME  FOR  SECOND  AIRFOIL  MOTION  TO  BEGIN. 

(O  TO  BEGIN  MOTION  AT  THE  SAME  TIME) . 

TOL  :  TOLERANCE  CRIERION  FOR  CONVERGENCE  FOR  (Uw)k  and  (Vw)k. 

TADJ  :  FACTOR  BY  WHICH  DTS  WILL  BE  ADJUSTED. 

SCLA  :  STEADY  LIFT  COEFF.  FOR  THE  SINGLE  AIRFOIL  AT  THE  SPEC.  AOA. 

SCM  :  STEADY  MOMENT  COEFF.  FOR  THE  SINGLE  AIRFOIL. 

SGAM  !  STEADY  VORTICITY  STRENGTH  FOR  THE  SINGLE  AIRFOIL. 

NGIES  ;  OPTION  TO  CHANGE  THE  UNSTEADY  KUTTA  CONDITION. 

0  EQUAL  PRESSURE  AT  THE  TRAILING  EDGE  PANELS. 

1  EQUAL  TANGENTAL  VELOCITIES  AT  THE  TRAILING  EDGE  PANELS. 


Figure  3.2  USPOTF2F  Input  Namelist  for  Flutter  Analysis  (Modified) 
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Cl  Magnitude  (1.0  deg  pitch;  pivot 
about  LE;  NACA0007;  200  panels; 

3  cycles,  65  time  steps  per  cycle) 


Figure  3.3  Cl  Magnitude  for  Pitch  about  Leading  Edge 
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Cl  Phase  (1  deg  pitch;  pivot  about 
LE;  NACA0007;  200  panels;  3 
cycles,  65  time  steps  per  cycle) 


Figure  3.4  Cl  Phase  for  Pitch  about  Leading  Edge 


Cm  Magnitude 


Cm  Magnitude  (1.0  deg  pitch;  pivot 
about  LE;  NACA0007;  200  panels; 

3  cycles,  65  time  steps  per  cycle) 


Figure  3.5  Cm  Magnitude  for  Pitch  about  Leading  Edge 
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Cm  Phase  (rad) 


Cm  Phase  (1  deg  pitch;  pivot  about 
LE;  NACA0007;  200  panels;  3 
cycles,  65  time  steps  per  cycle) 


Figure  3.6  Cm  Phase  for  Pitch  about  Leading  Edge 
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TABLE  3. 1  PITCH  OSCILLATION  PIVOT  L.E. 


Cl  Magnitude 


Cl  Magnitude  (1  deg  pitch;  0.25c 
pivot;  NACA0007;  200  panels;  3 
cycles,  65  time  steps  per  cycle) 


Figure  3.7  G  Magnitude  for  Pitch  about  1/4  Chord 


Cl  Phase  (1  deg  pitch;  0.25  pivot; 
NACA0007;  200  panels;  3  cycles 
65  time  steps  per  cycle) 


Figure  3.8  Cl  Phase  for  Pitch  about  1/4  Chord 


Cm  Magnitude 


Cm  Magnitude  (1  deg  pitch;  .25c 
pivot;  NACA0007;  200  panels; 

3  cycles  65  time  steps  per  cycle) 


Figure  3.9  Cm  Magnitude  for  Pitch  about  1/4  Chord 
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Cl  Magnitude  (1  deg  pitch;  0.37c 
pivot;  NACA0007;  200  panels; 

3  cycles,  65  time  steps  per  cycle) 


1 


Figure  3.11  Cl  Magnitude  for  Pitch  about  0.37c 


Cl  Phase  (rad) 


Figure  3.12  Cl  Phase  for  Pitch  about  0.37c 


Cm  Magnitude 


Cm  Magnitude  (1  deg  pitch;  0.37c 
pivot;  NACA0007;  200  panels;  3 
cycles,  65  time  steps  per  cycle) 


Figure  3.13  Cm  Magnitude  for  Pitch  about  0.37c 


Cm  Phase  (rad) 


Cm  Phase  (1.0  deg  pitch;  0.37c 
pivot;  NACA0007;  200  panels;  3 
cycles,  65  time  steps  per  cycle) 

o 

-0.1 
-0.2 
-0.3 
-0.4 
-0.5 
-0.6 
-0.7 
-0.8 
-0.9 
-  1 

0  2  4  6  8 

Kp  =  2bw/U 

Figure  3.14  Cm  Phase  for  Pitch  about  0.37c 
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Cl  Magnitude 


Cl  Magnitude  (0.01  h/2b  plunge; 
moment  about  LE;  NACA0007; 
200  panels;3  cycles,  65  time 
steps  per  cycle) 


Figure  3.15  Cl  Magnitude  for  Plunge 


Cl  Phase  (rad) 


Cl  Phase  (.01  h/2b  plunge;  moment 
about  LE;  NACA0007;  200  panels; 

3  cycles,  65  time  steps  per  cycle) 


Figure  3.16  Cl  Phase  for  Plunge 
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Cm  Magnitude 


Cm  Magnitude  (0.01  h/2b  plunge; 
moment  about  LE;  NACA0007;  200 
panels;3  cycles,  65  time  steps 
per  cycle) 


0  2  4  6  8 


Kp  =  2bw/U 


Figure  3.17  Cm  Magnitude  for  Plunge 


Cm  Phase(.01  h/2b  plunge;  moment 
about  LE;  NACA0007;  200  panels; 

3  cycles,  65  time  steps  per  cycle) 


Figure  3.1 8  Cm  Phase  for  Plunge 


IV.  SINGLE  DEGREE  OF  FREEDOM  FLUTTER 

When  discussing  flutter  of  a  typical  airfoil  section,  many  text  books  begin  the 
analysis  of  the  flutter  problem  as  a  two-dimensional  flexure-torsion  problem. 
Smilg  [Ref.  6]  presents  a  theoretical  study  for  one  degree  of  freedom  instability  of 
pitching  oscillations  of  an  airfoil  in  incompressible  flow.  Theoretical  results  show 
that  the  single  degree  of  freedom  instability  may  occur  if  the  airfoil’s  moment  of 
inertia  is  sufficiently  large  and  if  the  pitch  axis  is  located  forward  of  the  quarter 
chord.  The  pitch  axis  should  not  be  too  far  forward  of  the  leading  edge  so  that 
the  wing  becomes  stable.  It  is  important  to  distinguish  between  stall  flutter  due  to 
flow  separation  and  inviscid  incompressible  attached  flow  flutter  as  discussed 
here. 

A.  EQUATION  OF  MOTION  FOR  STEADY-STATE  OSCILLATIONS  IN 
ONE  DEGREE  OF  FREEDOM. 

The  classical  spring-mass  damper  equation  of  motion  is: 

mx  +  cx  +  kx  =  F 

The  analogous  equation  for  the  pitch  degree  of  freedom  of  a  rigid  airfoil  in  a  two 
dimensional  flow  is: 

Ia  +  ga+  ka  =  Ma 

Assuming  no  structural  damping  (g=0),  k  =  I©2  ,  a  =  -(o2a0e,a>\  and 
a  =  a0eia* : 

I(co2-a>2)  =  -Ma 

Defining  the  aerodynamic  pitching  moment  as  Ma  =  Cma^-pU24b2  and 
substituting  into  the  above  equation  gives: 
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For  undamped  oscillations  of  the  airfoil  to  exist,  the  non-dimensional  inertia 
parameter  I0/7tpb4  must  exceed  approximately  550.  As  the  elastic  stiffness 

increases,  the  inertia  parameter  must  also  increase  for  instability  to  occur.  For  a 
full  scale  aircraft  at  sea  level,  the  inertia  parameter  typically  will  not  exceed  30  for 
control  surfaces,  wings  or  stabilizers  [Ref.  6]. 

For  the  oscillations  to  damp  out,  Imaginary(CMa)  must  be  negative  and  the 

net  work  done  on  the  airfoil  per  cycle  must  be  negative.  Negative  work  implies 
that  the  airfoil  is  performing  work  on  the  air,  dissipating  energy  from  the  wing. 
When  the  aerodynamic  moment  and  pitch  displacement  are  in  the  same  direction 
throughout  the  entire  cycle  positive  work  is  performed  on  the  airfoil.  This 
transfer  of  energy  is  the  origin  of  negative  damping  which  leads  to  instability  of 
pitch  oscillations.  The  neutral  stability  boundary  corresponds  to  180  degrees 
phase  lag  between  airfoil  motion  and  lift  and  moment  time  history.  This  phase  lag 
is  defined  by  the  necessary  condition: 

Imag(CMo)  =  0 

The  real  part  of  the  lift  and  moment  coefficients  represents  the  in-phase 
aerodynamic  loading  while  the  imaginary  part  of  the  coefficients  is  the  out-of- 
phase  portion. 
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B.  SINGLE  AIRFOIL  FLUTTER  VERIFICATION 

Pitch  damping  (Imag(CMo))  is  plotted  for  several  amplitudes  in  Figure  4.2. 

The  critical  reduced  frequency  representing  the  flutter  boundary  is  approximately 
0.115  compared  to  a  theoretical  value  of  0.078.  The  accuracy  of  the  critical 
reduced  frequency  increases  as  the  pitch  amplitude  is  decreased. 

C.  AIRFOIL  IN-GROUND  EFFECT  SIMULATION 

To  perform  a  ground  effect  simulation,  the  additional  boundary  condition  of 
flow  tangency  at  ground  level  requires  a  second  airfoil’s  reflection  about  ground 
level.  Distance  from  the  ground  is  non-dimensionalized  by  chord  length.  First, 
the  steady  state  effects  were  examined  by  placing  the  two  airfoils  at  zero  degrees 
angle  of  attack  and  varying  the  vertical  distance  between  them.  At  less  than  3 
chord  lengths  of  separation,  the  venturi  effect  or  accelerated  flow  between  the 
airfoils  causes  a  decrease  in  lift  corresponding  to  a  positive  pitching  moment 
[Figure  4.3]. 

The  critical  reduced  frequency  of  a  pitching  airfoil  in  ground  effect  initially 
increases  slightly  then  decreases  rapidly  as  the  airfoil  separation  is  reduced 
[Figure  4.4],  For  an  airfoil  with  a  one  foot  chord,  a  natural  frequency  of  one  cycle 
per  second  and  inertia  parameter  i/rcpb4  equal  to  550,  the  panel  code  calculates 

the  flutter  boundary  to  occur  at  a  reduced  frequency  of  Kp  equal  to  0.115  and 
CMo  equals  -0.026613  out-of-ground  effect.  The  out-of-ground  effect  flutter 
speed  is  80  ft/sec  and  flutter  frequency  co/coa  equals  1.464.  At  4.5  chord  lengths 

separation  or  2.25  chord  lengths  above  ground  level  the  critical  reduced 
frequency  Kp  equals  0.101  and  CMa  equals  -0.028942,  which  corresponds  to  an 
in-ground  effect  flutter  speed  of  125  ft/sec  and  flutter  frequency  co/coa  of  2.011. 
Clearly,  the  flutter  speed  and  flutter  frequency  increase  in-ground  effect. 
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The  same  type  of  simulation  was  conducted  for  t\  o  airfoils  pitching  in-phase. 
The  flutter  speed  decreased  and  the  critical  reduced  frequency  of  the  flutter 
boundary  increased  [Figure  4.5]. 

A  physical  understanding  of  these  results  requires  visualization  of  the  wake. 
The  in  ground  effect  simulation  tends  to  cancel  the  wake’s  effect  because  the 
reflected  wake  has  the  same  strength  and  opposite  direction  [Figure  4.6],  In- 
phase  oscillations  will  enhance  the  wake’s  effect  upon  the  aerodynamics  of  the 
airfoil  [Figure  4.7].  The  two  wakes  maintain  the  same  separation  as  the  airfoils 
have  and  the  vortices  have  the  same  direction. 
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Figure  4.1  Phase  Diagram  for  Airfoil  Stability 
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Cm  Imaginary 


Im(Cm)  for  Various  Pitch 
Magnitudes  (pivot  about  LE; 
NACA0007;  200  panels;  3 
cycles,  65  time  steps  per  cycle) 


Figure  4.2  Single  Airfoil  Flutter  Verification 
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Steady  State  Flow  of  NACA  0007 
Airfoil  in  Ground  Effect;  0  deg 
AOA;  Moment  about  LE;  200  panels 


Figure  4.3  Steady  State  Aerodynamics,  Moment  Referenced  to  Leading  Edge 
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Figure  4.4  Flutter  Boundary  In  Ground  Effect 
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TABLE  4.1  REFERENCES  FIGURE  4.2 


Cm  Imaainary  due  to  Ditch 


NACA0007:  200  panels;  3  cycles  65  steps  per  cycle 


0.02 

0.000537 

0.001857 

0.003895 

0.04 

0.000639 

0.002247 

0.005146 

0.06 

0.000570 

0.002043 

0.004750 

0.08 

0.000393 

0.001456 

0.003459 

0.10 

0.000126 

0.000594 

0.001610 

0.12 

-0.000217 

-0.000480 

-0.000638 

0.14 

•0.000618 

-0.001718 

-0.003193 

0.16 

-0.001065 

-0.003084 

-0.005978 

0.18 

-0.001545 

-0.004544 

-0.008884 

0.20 

•0.002060 

-0.006101 

-0.012034 

TABLE  4.2  REFERENCES  FIGURE  4.3 
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TABLE  4.3  REFERENCES  FIGURE  4.4 


TABLE  4.4  REFERENCES  FIGURE  4.5 
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ase  Oscillations 


Figure  4.6  Out-of-Phase  Oscillations,  Wake  Cancellation 
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Dimensional  Time 
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Non-Dimensional  Time 


V.  ACTIVE  FLUTTER  CONTROL 


The  traditionally  accepted  methods  of  moving  an  airfoil’s  flutter  boundary  is 
shifting  the  center  of  gravity  or  stiffening  the  structure.  These  two  methods 
affect  the  inertial  and  elastic  forces.  In  the  past,  aerodynamic  forces  were 
accepted  as  an  uncontrollable  quantity  dependent  upon  reduced  frequency. 
Through  active  flutter  control,  these  forces  can  be  controlled. 

An  investigation  of  how  the  one-dimensional  flutter  boundary  may  be  shifted 
by  positioning  an  oscillating  control  airfoil  upstream  of  a  neutrally  stable 
reference  airfoil  is  presented  in  this  chapter.  The  aerodynamic  forces  may  be 
modified  to  stabilize  an  airfoil. 

A.  INVESTIGATION  OF  AIRFOIL  POSITIONING  AND  SIZE  FOR 

ACTIVE  FLUTTER  CONTROL. 

The  USPOTF2F  code  can  accommodate  and  scale  a  second  control  airfoil 
while  the  reference  airfoil  remains  fixed  at  the  origin  of  the  global  coordinate 
system.  The  variables  manipulated  are  horizontal  spacing,  vertical  spacing  and 
scaling  of  the  contro  foil.  The  physical  pitching  frequencies  of  the  two  airfoils 
are  equal  for  all  cast  presented  in  this  chapter.  The  pitch  axis  is  located  at  the 
leading  edge. 

Pitch  damping  results  are  dependent  upon  the  wavelength  of  the  wake, 
defined  as: 

v  _2nc 
Kp  X 

where  Kp  is  reduced  frequency  based  on  full  chord,  c  is  chord  length  and  X  is 
the  wavelength  of  the  wake.  Pitch  damping  refers  to  the  imaginary  part  of  the 


66 


moment  coefficient.  Positive  Imag(Cm)  corresponds  to  negative  pitch  damping  or 
instability. 

Figure  5.1  displays  the  pitch  damping-frequency  dependence  for  identically 
sized  airfoils  staggered  by  two  chord  lengths  ahead  and  below  the  reference 
airfoil.  Figure  5.1  is  very  similar  to  Figure  4.2  constructed  for  single  airfoil  flutter 
verification.  The  staggered  airfoil  curve  [Fig.  5.1]  shifts  down  and  slightly  to  the 
left  when  compared  to  Figure  4.2. 

Figure  5.2  shows  the  airfoil  interaction  at  positive  and  negative  Xshift  values. 
For  positive  Xshift,  the  pitch  damping  asymptotically  approaches  the  single  airfoil 
value  of  0.000126  [Table  4.1].  Between  Xshift  values  of  -5  to  +10  the  airfoil 
influence  supersedes  wake  effects.  The  periodicity  of  pitch  damping-Xshift  is  the 
wake  wavelength.  The  mean  of  the  sine  curves  in  Figure  5.3  through  Figure  5.5 
follows  the  single  airfoil  pitch  damping  curve  [Fig.  4.2].  As  the  reduced 
frequency  increases,  the  amplitude  of  Imag(Cm)  decreases  and  the  airfoil  becomes 
more  stable.  In  summary,  the  frequency  dependence  of  pitch  damping  is  inversely 
proportional  to  Cm  imaginary  period,  mean  and  amplitude. 

The  effect  of  scaling  or  sizing  the  control  airfoil  is  displayed  in  Figures  5.6 
and  5.7.  The  ability  of  the  control  airfoil  to  affect  the  flutter  boundary  decreases 
proportionally  with  control  airfoil  size.  As  the  size  of  an  airfoil  decreases,  the 
strength  of  the  wake  core  vortex  shed  from  the  trailing  edge  also  decreases  due 
reduced  circulation  around  the  control  airfoil. 

B.  APPLICATION  TO  UNSTEADY  AERODYNAMICS  OF  ROTARY 
WINGS. 

The  USPOTF2F  code  is  well  suited  for  unsteady  aerodynamic  simulation  of  a 
hovering  helicopter.  A  two  dimensional  approximation  considers  the  influence  of 
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shed  vortices  that  have  translated  below  the  rotor  disc  and  passed  under  the  next 
blade.  Wake-excited  flutter  has  been  documented  in  laboratory  tests  and  during 
initial  whirl  tower  testing  of  the  AH-56  Cheyenne.  Wake-excited  flutter  will  only 
occur  at  low  inflow  rates  when  wake  spacing  is  minimal.  The  low  inflow 
configuration  occurs  during  ground  operations  when  blade  pitch  angles  are  small. 
In  an  actual  helicopter  rotor  flow  field,  the  shed  vorticity  sheets  will  superimpose 
upon  each  other  becoming  stronger  with  the  passage  of  each  successive  blade. 
A  drawback  to  the  USPOTF2F  code  is  that  it  only  accounts  for  the  wake  of  one 
previous  blade. 

C.  COMPARISON  TO  TWO  DIMENSIONAL  APPROXIMATION  TO 

WAKE  FLUTTER  OF  ROTARY  WINGS  BY  ROBERT  LOEWY. 

1.  Parameters 

In  Reference  7,  page  90,  Figure  15  Loewy  presents  a  graph  of  pitch 
damping  coefficient  versus  frequency  ratio  as  a  function  of  inflow  parameter  [Fig. 
5.8].  The  inflow  parameter  is  the  vertical  spacing  of  the  wakes  of  previous  blades 
corresponding  to  Yshift  in  USPOTF2F  code.  The  frequency  ratio  m=co/Q  is 
divided  into  two  portions,  the  integer  and  non-integer  part. 

The  non-integer  portion  of  m  represents  wake  phasing.  When  m  equals 
0.5,  this  corresponds  to  the  wakes  being  out  of  phase  as  shown  in  Figure  4.6. 
The  vortex  from  the  control  airfoil  or  preceding  blade  is  directly  below  a  core 
vortex  from  the  reference  airfoil  in  the  opposite  direction.  The  wakes  are  counter 
phase  or  1 80  degrees  out  of  phase.  When  m  equals  zero  or  one,  the  phase  shift  is 
0  and  360  degrees  respectively  [Fig.  4.7].  The  vortex  wake  sheet  has  equal 
spacing  vertically  with  equal  vortex  strength  and  direction  for  identical  airfoils. 
Couch  in  Reference  10,  page  39  treats  m  exclusively  as  a  phasing  parameter  in  his 
finite  wake  theory. 
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The  integer  portion  of  m  represents  the  reduced  frequency  of  the  airfoil. 


o) 


Multiplying  by  the  inverse  of  aspect  ratio  for  a  rectangular  wing  yields  reduced 
frequency. 


0)  c  cue  c  c 


Reduced  frequency  is  physically  defined  as  the  ratio  of  chord  length  to  wake 
wavelength  as  shown  above. 

2.  Comparison  to  Loewy  Results 

Figure  5.8  shows  areas  of  negative  pitch  damping  between  m=0.5  to  1.0 
and  m=1.5  to  2.0.  The  m=0.5  to  1.0  region  corresponding  to  a  phase  shift  of  180 
to  360  degrees  is  illustrated  in  Figure  5.3  when  Xshift  ranges  from  -31.440  -62.8, 
which  is  a  region  of  instability.  Figure  5.3  through  Figure  5.5  confirms  that  as  the 
reduced  frequency  increases  (Kp=0.1  to  1.0),  the  amplitude  of  the  pitch  damping 
curve  will  decrease.  The  difference  in  shape  of  Loewy’s  pitch  damping  curve 
[Fig.  5.8]  and  sinusoidal  shape  of  Figure  5.3  is  due  to  Loewy’s  use  of  an  infinite 
number  of  wakes  in  his  classical  closed  form  solution  while  USPOTF2F  code  only 
accounts  for  one  wake. 

D.  USPOTF2F  CODE  LIMITATIONS 

Theodorsen’s  flutter  theory  for  simple  harmonic  motion  assumes  sinusoidally 
varying  lift  and  moment  coefficients,  which  are  represented  mathematically  by 
e,u* .  The  primary  limitation  is  that  the  unequal  physical  frequencies  of  control 
and  reference  airfoil  yield  higher  harmonic  time  histories  of  sine  and  cosine.  The 
phase  subroutine  attempts  to  curve  fit  a  sine  wave  to  a  higher  harmonic  time 
history  resulting  in  rapid  inconsistent  phase  shifts. 
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Cm  Imaginary 


Cm  Imaginary  Values  for 
Staggered  Airfoils 
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Figure  5.1  Frequency  Dependence  of  Pitch  Damping  Xshift=-2,  Yshift=-2 
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TABLE  5.1  DATA  REFERENCES  FIGURE  5.1 


Pitching  Airfoils  in  Phase  at  Various  Freauencies 


Pivot  about  LE:  NACA0007;  200  panels;  3  cycles  65  ste 


Xshift  -2.0;  Yshift  -2.0;  Oaipl  -  1.0;  Dalp2  -  1.0; 


Airfoils  1  and  2  freauencv  sweep  from  .025  to  .250; 


Freal 


0.025 


0.038 


0.050 


0.063 


0.075 


0.100 


0.113 


0.125 


0.138 


0.150 


0.163 


0.175 


0.188 


0.200 


0.213 


0.225 


0.238 


0.250 


Frea2 


0.025 


0.038 


0.050 


0.063 


0.075 


0.088 


0.100 


0.113 


0.125 


0.138 


0.150 


0.163 


0.175 


0.188 


0.200 


0.213 


0.225 


0.238 


0.250 


m 


Cm  Imaaina 


0.000959 


0.000972 


0.000822 


0.000545 


0.000222 


SFTimn 


-0.000549 


.000982 


-0.001432 


-0.001887 


-0.002345 


•0.00282 


.003297 


-0.003774 


•0.004234 


-0.004698 


•0.005149 


•0.005603 


-0.00607 
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Imaginary  Cm 


Cm  Imaginary  Plot  for  Constant  Kp 

0.1 

0.002 
0.0015 
0.001 
0.0005 
0 

-0.0005 
-0.001 
-0.0015 
-0.002 
-0.0025 
-0.003 

-  1  0  0  1  0  20 
Xshift  (chord  lengths) 

Figure  5.2  Pitch  Damping  at  Small  Xshift  Values  Compared  to  Wake  Wavelength 
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TABLE  5.2  DATA  REFERENCES  FIGURE  5.2 


Cm  Imaginary  Plot  for  Constant  Kp 


Figure  5.3  Pitch  Damping-Xshift  Dependence 


Imaginary  Cm 


Cm  Imaginary  for  Constant  Kp  =  0.2 


Figure  5.4  Pitch  Damping-Xshift  Dependence 
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TABLE  5.3  DATA  REFERENCES  FIGURE  5.3  AND  FIGURE  5.4 


Pitching  Airfoils 


Pivot  about  LE;  NACA0007:  200  panels;  3  cycles  65  steps  per 


Yshift  -2.0:  Dalpl  1.0;  Dalp2  1.0 


Xshift  varied  from  0.0  to  -120  chord  lengths 


-  0.1 


Xshift 


-6.28 


-12.57 


-18.85 


-25.13 


-31.42 


-37.70 


-43.98 


-50.26 


-56.55 


-62.83 


-69.11 


-75.40 


-81.68 


-87.96 


-94.25 


100.53 


106.81 


113.10 


119.38 


Cm  Imagina 


0.000737 


-0.001992 


•0.002929 


•0.002595 


-0.001021 


0.001010 


0.002836 


0.003599 


0.003089 


0.001560 


-0.000474 


-0.002236 


-0.003109 


-0.002814 


-0.001184 


C. 000944 


0.002751 


0.003464 


0.002851 


0.001307 


ni 


K 


Xshift 


0.00 


-3.14 


-6.28 


-9.42 


-12.57 


-15.71 


-18.85 


-21.99 


-25.13 


-28.27 


-31.42 


-34.56 


-37.70 


-40.84 


-43.98 


-47.12 


-50.27 


-53.41 


-56.55 


-59.69 


Cm  Imagina 


-0.001999 


-0.005072 


-0.006196 


-0.005637 


-0.003424 


-0.000595 


0.002647 


0.001769 


-0.000400 


-0.003206 


-0.005419 


-0.006422 


-0.005997 


-0.003724 


-0.000667 


0.001671 


0.002466 


0.001482 


-0.00061 1 
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TABLE  5.4  DATA  REFERENCES  FIGURE  5.5 


Pitchina  Airfoils  Ko  -  t.O 


Pivot  about  LE;  NACA0007;  200  panels; 


Yshift  *2.0;  Dalpl  1.0;  Dalp2  1.0 


Xshift  varied  from  0.0  to  -18.22  chord  lengths 


Xshift 


0.00 


-0.63 


-1 .26 


-1 .88 


-2.51 


-3.14 


-3.77 


-4.40 


-5.03 


-5.65 


-6.28 


-6.91 


-7.54 


-8.17 


•8.80 


-9.42 


-10.05 


-10.68 


-11.31 


-11.94 


-12.57 


-13.19 


-13.82 


-14.45 


-15.08 


-15.71 


-16.34 


-16.96 


-17.59 


-18.22 


5  cycles  39  ste 


-0.027219 


-0.0280  60 


-0.029022 


-0.029452 


.029082 


-0.027958 


•0.026568 


-0.025448 


-0.025049 


-0.025478 


•0.026566 


-0.027912 


-0.029073 


-0.029587 


-0.029219 


wznm 


-0.026656 


-0.025488 


•0.025106 


-0.025538 


-0.026627 


-0.027934 


-0.029100 


-0.029640 


-0.029299 


-0.028121 


.026633 


-0.025452 


.025088 


-0.025537 


3  cycles  65  ste 


-0.027008 


-0.027836 


-0.028742 


•0.029159 


•0.028802 


-0.027730 


-0.026411 


-0.025333 


-0.024928 


-0.025330 


-0.026357 


-0.027685 


-0.028777 


-0.029245 


•0.028880 


-0.027792 


mmn 


-0.025244 


-0.024793 


-0.025154 


m 


TABLE  5.5  DATA  REFERENCES  FIGURE  5.6 


Pitching  Airfoils  Kp  -  0.1 


Pivot  about  LE;  NACA0007;  200  panels:  3  cycles  65  steps  per  cycle; 


Yshift  -2.0:  Daipl  =  1.0:  0alp2  «  1.0: 


Xshift  varied  from  0  to  120  chord  lengths  for  scales  of  0.25.  0.50  and  1 


ml 


Xshift 


0.00 


•6.28 


-12.57 


-18.85 


-25.13 


-31.42 


7.70 


-43.98 


-50.26 


-56.55 


-62.83 


-69.11 


-75.40 


-81.68 


-87.96 


-94.25 


100.53 


106.81 


119.38 


Scale  1 .0 


0.000737 


•0.001992 


-0.002929 


-0.002595 


-0.001021 


0.001010 


0.002836 


0.003599 


0.003089 


0.001560 


.000474 


-0.002236 


-0.003123 


-0.002814 


-0.001184 


0.000944 


0.002751 


0.003464 


0.002851 


0.001307 


Scale  0.5 

Scale  0.25 

0.000429 

0.000317 

-0.000939 

-0.000383 

-0.001471 

-0.000638 

-0.00125 

-0.000532 

-0.000415 

-0.000114 

0.000659 

0.000437 

0.000903 

0.002005 

0.001116 

0.001742 

0.000988 

0.000901 

0.000557 

-0.000172 

0.000013 

•0.001098 

-0.000455 

-0.001599 

-0.000697 

-0.001355 

-0.000581 

-0.000495 

-0.000154 

0.000597 

0.000402 

0.001536 

0.000869 

0.001932 

0.001076 

0.001641 

0.000929 

0.000765 

0.000477 
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Cm  Imaginary 


Yshift  Dependence  Upon  Pitch 
Damping  (Xshiffe-43.98) 


Figure  5.7  Yshift  Dependence 
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TABLE  5.6  DATA  REFERENCES  FIGURE  5.7 


Figure  5.8  Loewy  Pitch  Damping  Curves  [Ref.  7] 
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VI.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  SINGLE  AIRFOIL  ANALYSIS 

Aerodynamic  verification  of  UPOTFLUT  code  is  consistent  with  results  by 
Riester  [Ref.  1].  Incompressible  bending-torsion  flutter  computations  compare 
favorably  with  Theodorsen  theory  when  the  natural  frequency  ratios  are  less  than 
1.2.  It  was  shown  that  the  critical  reduced  frequency  may  differ  by  up  to  a  factor 
of  three  from  Theodorsen’s  result  for  high  critical  frequencies  and  high  natural 
frequency  ratios. 

B.  TWO  AIRFOIL  ANALYSIS 

The  nonlinear  theory  for  simple  harmonic  motion,  magnitude  and  phase 
relationships  that  exist  between  airfoil  motion  and  aerodynamic  forces  have  been 
verified  by  comparison  to  the  classical  Theodorsen  analysis.  Agreement  with 
theory  consistent  with  previous  panel  codes  was  obtained.  Aerodynamic 
coefficient  magnitude  errors  were  less  than  five  percent  and  phase  errors  typically 
ranged  from  one  to  three  degrees. 

A  single  degree  of  freedom  flutter  analysis  of  an  airfoil  in-ground  effect 
simulation  showed  that  the  flutter  speed  increases  as  ground  distance  decreases 
due  to  wake  cancellation.  The  one  dimensional  flutter  boundary  may  be  shifted 
by  positioning  and  oscillating  a  control  airfoil  upstream  of  a  neutrally  stable 
reference  airfoil.  An  investigation  of  positioning  and  size  of  a  control  airfoil  for 
active  flutter  control  showed  the  pitch  damping  was  periodic  with  wake  phasing. 
Pitch  damping  period,  mean  and  amplitude  are  inversely  proportional  to  reduced 
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frequency.  A  trend  comparison  with  the  wake  flutter  theory  by  Loewy  illustrated 
consistent  code  results. 

Theodorsen  flutter  theory  for  simple  harmonic  motion  assumed  sinusoidally 
varying  lift  and  moment  coefficients,  which  were  represented  mathematically  by 
e10*.  The  physical  frequency  differences  between  the  control  and  reference 
airfoil  yielded  higher  harmonics  in  the  time  history,  causing  Theodorsen  theory  to 
be  invalid. 

C.  RECOMMENDATIONS 

A  drawback  to  this  code  was  the  lengthy  computer  run  times  which  were 
approximately  60  hours  per  figure  on  a  100  Mhz  SGI  Workstation.  Run  times 
could  be  reduced  by  running  the  code  in  parallel.  Dividing  the  simulation  into 
separate  directories  or  cases  and  running  up  to  ten  workstations  simultaneously 
make  time  requirements  manageable.  This  requires  the  availability  of  batch 
queues  and  access  to  a  network.  Streamlining  of  the  code  by  an  experienced 
programmer  is  recommended. 

The  frequency  dependence  of  the  control  airfoil  upon  a  reference  airfoil  was 
not  resolved.  The  non-sinusoidal  lift  and  moment  time  histo.  .  will  require 
computations  to  be  completed  in  the  time  domain.  The  concept  of  active  flutter 
control  is  feasible  and  certainly  worthy  of  further  theoretical  study. 
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APPENDIX 


This  appendix  contains  theoretical  Theodorsen  flutter  analysis  as  described 
by  Fung  [Ref.  8].  The  MATLAB  program  was  written  for  the  Macintosh.  Next, 
time  histories  of  lift,  moment  and  pitch  angle  for  various  cases  are  illustrated. 
Finally,  three  sample  wake  position  plots  show  the  core  vortex  position  at  the 
final  time  step. 
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'  11us  Program  performs  theoretical  theocicirsieh  analysis  as 

%  described  hy  Fir^. 

clear 

clg 

% 

%  Aeroc^namic  Coefficients 
K(l)  =  100; 

Ui(l)  =  1; 

Maipha(l)  =  0.375; 

Lalpha(l)  =  0.5; 

K(2)  =  4; 

LJi(2)  =  0.9848-0.2519*1; 

Maipiia(2)  =  0. 375-0. 250*i; 
tnlpha(2)  =  0.42179-0.49423*1; 

K (3 )  =  2; 

Lh(3)  =  0.9423-0.5129*1; 
tolplia(3)  =  0.375-0.5*1; 

Lalplia(3)  =  0.18580-0.98405*1; 

K(4)  =  1.2; 

Ui(4)  =  0. 8538-0. 8833*i; 

Malpha(4)  =  0.375-0.83333*1; 

Lalpha(4)  =  -0.38230-1.59487*1; 

K(5)  =  0.8; 

Ui(5)  =  0.7088-1.3853*1; 

*felpha(5)  s  0.375-1.25*1; 

Lalpha(5)  =  -1.5228-2.27119*1; 

K(6)  =  0.6; 

Ui(6)  =  0.5407-1.9293*1; 
telpha(6)  =  0.375-1.66667*1; 
talpta(6)  =  -3.17490-2.83045*1; 

K(7)  =  0.5; 

Ui(7)  =  0.3972-2.3916*1; 

*felpha(7)  =  0.375-2.0*1; 

Lalpha(7)  =  -4.8860-3.1860*1; 

K(8)  =  0.4; 

Ui(Q)  =  0.1752-3.1250*1; 
ffelrJa(8)  =  0.375-2.5*1; 

Lalpha(8)  =  -8.1375-3.5625*1; 

K(9)  =  0.34; 

Ui(9)  =  -0.022-3.8053*1; 
tolplia(9)  =  0.375-2.94118*1; 

Lalpha(9)  =  -11.714-3.7396*1; 

K(10)  =  0.3; 

th(10)  =  -0.1950-4.4333*1; 
tolpiH(lO)  =  0.375-3.333333*1; 
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tnlpha(lO)  =  -15. 473-3. 7822*i; 

K(ll)  =  0.266666667; 

Lh(ll)  =  -0.3798-5. 1084*i; 
tolpha(ll)  =  0. 375-3. 75*i; 

Lnlpha(ll)  =  -20. 0337-3. 6847*i; 

K(12)  =  0.24; 

Iii(12)  =  -0. 552-5. 8242*i; 

Nftlp)n(12)  =  0.375-4.16667*1; 

Lalplia(12)  =  -25. 319-3. 5256*i; 

K(13)  =  0.2; 

Hi  (13)  =  -0. 886-7. 2760*i; 

Maltfia(i3)  =  0. 375-5. 0*i; 

Lalplra(13)  =  -37.766-2.8460*1; 

Mi  =  0.5; 

%  Ncn-dimaisional  coefficients 

q  =  0; 
vlwa  =  1.6; 
mi  =  4; 
a  =  -0.3; 
xaliiia  =  0.2; 
ralfiia2  =  0.25; 

%  Flutter  Determinant 

%  Search  for  real  root  using  bisection  method 
%  Outside  for  loop  conducts  the  frequency  sweep 


for  m  =  1:13 

X(l)  a  0.2; 

X(3)  =  2; 

X (2)  =  (X(3)-X(l) )/2; 
for  its  =  1:50 

A  =  nit*(l-X*(vAiwa~2)*(Ui*g))dii(m); 

B  =  un*xalplia>Lalpha(m)-Lli(m)*(0.5+a); 

D  =  nu*xalpha+Mi-Lh(m)*(0.5+a) ; 

E=iiii*ralf)ha2*  (l-X*(l+i*g)  )-0.5*(0.5+a)4felpha(m)-Laipha(m)  *(0.5+a)+Iii(m)  *(0.5+a 

"2; 

DET  =  A.  *E  -  B*D; 

R  a  real  (DET)  ; 
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if  R(1)*R(2)  <  0 
X{3)  =  X(2) } 
else 

X(l)  =  X(2) } 
aid 

X(2)  =  (X(3)4X(l))/2; 
ati 

Rroot  (m)  =  X(2); 

Rerror(m)  =  X(3)  -  X(l)‘ 

%  Search  for  imaginary  root 

X(l)  =  0.2; 

X(3)  =  2; 

X(2)  =  (X(3)~X(1)  )/2{ 
for  its  =  1:50 


A  =  mi*  (1-x*  (whwa/v2)  *  (l4i*g) )  +Lh(m) ; 

B  =  nu*xalpha+Lalpha(m)  -Lh(m) *  (0.5+a)  * 

D  s  rru*xalpha+fti-Lh (mj  *  (0 . 54a) ; 

E=*u*ralpha2*  (1-x*  (Ui*g) )  -0 . 5*  (0 . 54a)  4Maipha  (in)  -Ulfcha  (m)  *(0  <  5+a)  4lii  (m)  *  (0 . 54a 
~2; 

DBT  =  A.*E  -  B*D; 

I  =  inBg(DET)  ; 
if  1(1) *1(2)  <  0 
X(3)  =  X(2)  ,* 
else 


X(l)  =  X(2) ; 
ax! 

X(2)  =  (X(3)4X(1) )/2) 
aid 

Iroot  (m)  =  X{2),‘ 

Ierror(m)  =  x(3)  -  X(l)  ,* 

end 

Rtrot 

Iroot 

cliff  =  Rroot  -  Iroot; 

%  Plot  of  real  and  imaginary  roots 


4 


%plot(l  ./  K,  Rroot, 1  ./  K,  Iroot) 
%grid 

fcdabelCl/k1) 
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%y label ( ' (X) * ) 

"tpause 

%  Search  of  roots  for  critical  frequency 


n  =  length  (Rroot),‘ 
for  m=l:n 

if  diff(m)  >  0 
krf=K(m) ; 
break 
end 
ad 

%  Calculation  of  the  flutter  ooefficent  that  lies  between 
*  m  and  m-1  reduced  frequency. 

IovK(l)  =  1/K(m-1); 

Xreal(l)  =  sqrt  (Rroot  (m-1) )  ,* 

Xiirag(l)  =  sqrt (Iroot (m-1) ) } 

IovK(2)  =  1/K(m)  ; 

Xreal(2)  =  sqrt  (Rroot  (m))‘ 

Xirag(2)  =  sqrt  (Iroot  (m) )  ‘ 

xl  =  IcwK(l)  :0.0001‘IovR(2)  ,• 
n  =  length  (xl) ; 

yl  =  l/(IovK(2)  -  lcwK(l) )  *  ( (IovK(2)  -  xl)*Xredi(l)  +  (xl  -  Iovfc(l))*Xreal(2))  ‘ 
y2  =  l/(IovK(2)  -  IovK(l) )*{(IovK(2)  -  xl) *Ximag(l>  +  (xl  -  Iovk(l))*XijiHg(2)); 

plot  (xl,yl,xl,y2) 
grid 

xlabel ( ' 1/k' ) 
y label ('sqrt(X) ') 


for  m  =  l:n 

if  y2(m)-yl(m)  <  0, 
Iok  =  xl(m) 
Xroot  =  y2  (m) 
break 
ad 
ad 

kcrit  *  1/Iok 
coef  =  Iok  /  Xroot 


* 
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about  LE;  Yshift=50 
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Nondimensional  Time 


Nondimensional  Time 


Nondimensional  Time 


Nondimensional  Time 
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Nondimensional  Time 


Scale  =  0.5 


Vortex  Position 


Vortex 


Wake  Trace  Kp(l)  =  1.0  Kp(2)  =  2.0  Scale  =  0.5 


Vortex  Position 
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